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ON THE DISCRETIZATION OF THE 
EULER-POINCARE-SUSLOV EQUATIONS IN 50(3) 


FERNANDO JIMENEZ AND JURGEN SCHEURLE 


Abstract. In this paper we explore the discretization of Euler-Poincare-Suslov 
equations in 50 (3). We prove that the consistency order of the unreduced and 
reduced setups, when the discrete reconstruction equation is given by a Cayley 
retraction map, are related to each other in a nontrivial way. Moreover, we 
give precise conditions under which general and variational integrators gen¬ 
erate a discrete flow preserving the distribution. These results are carefully 
illustrated by the example of the Suslov problem in 50(3), establishing general 
consistency bounds and illustrating the performance of several discretizations 
through some plots. Finally, we show that any constraints-preserving dis¬ 
cretization may be understood as being generated by the exact evolution map 
of a time-periodic non-autonomous perturbation of the original continuous¬ 
time nonholonomic system. 


1. Introduction 

The Lagrangian formulation of mechanical systems with nonholonomic cons¬ 
traints has been extensively studied in recent years (see [Sj HH [26] for a complete 
description and extensive bibliographies). In short, this kind of systems are cha¬ 
racterized by so-called nonholonomic constraints, i.e. constraints involving both 
configuration as well as velocity variables, and which can not be integrated to 
purely configuration-dependent constraints (in this case the constraints are called 
holonomic). Moreover, the preservation of the nonholonomic constraints by suitable 
integrators is a delicate issue upon which a lot of attention has been put by the Geo¬ 
metric Mechanics comunity (see the recent works, such as [ini miiiiiia i2ii2iiMj , 
which have introduced numerical integrators for nonholonomic systems with very 
good energy behavior and preservation properties). The approaches in most of these 
references are based on the ideas of |33l|38]. In these works, the continuous varia¬ 
tional principles are replaced by discrete ones aiming to obtain proper integrators 
approximating the continuous dynamics; we will call the integrators related to this 
framework variational integrators. Analogously, in the case of nonholonomic me¬ 
chanics, the continuous Lagrange-d’Alembert principle, which provides the actual 
dynamics, is replaced by a discrete Lagrange-d’Alembert principle on the discrete 
phase space. Of special interest are the seminal works on nonholonomic integration 
[ 13131 ], where a discrete version of the Lagrange-d’Alembert principle is proposed 
by introducing a proper discretization of the nonholonomic distribution. 

Reduction theory is one of the fundamental tools in the study of mechanical 
systems with symmetry and it concerns the removal of symmetries and the associ¬ 
ated conservation laws (see [32] for more details). Roughly speaking, we will say 
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that a Lagrangian system is invariant under a Lie group if the Lagrangian function 
is invariant under the tangent lift of the action of the Lie group on the config¬ 
uration manifold. Furthermore, for a symmetric mechanical system the process 
of reduction eliminates the directions along the group variables and thus provides 
a system with fewer degrees of freedom. When the configuration manifold is the 
Lie group itself, the case this work is concerned about, the process of reduction 
for unconstrained Lagrangian mechanics leads from the Euler-Lagrange equations 
to the Euler-Poincare ones. Needless to say, nonholonomic systems may also pos¬ 
sess symmetries, and the geometrical treatment of such situations was studied in 
[siiHiiniiii] among other references. On the other hand, variational integrators 
for reduced systems were carefully studied from the theoretical point of view in 
[28112^ . The combination of these two elements, say symmetric nonholonomic La¬ 
grangian systems and variational integrators, was studied in m for a ll system, 
i.e. a nonholonomic system whose configuration manifold is the Lie group G, and 
where both the Lagrangian and the constraint distribution are invariant with re¬ 
spect to the induced left action of G on TG. The dynamics of this kind of system 
is described by the Euler-Poincare-Suslov equations [23 EH Eg. Other interesting 
discretizations of symmetric nonholonomic systems can be found in m- 

The Euler-Poincare-Suslov equations on 50(3) are the object of study in this 
paper, more concretely from the point of view of the results presented in ^D] . 
As already pointed out, the Euler-Poincare-Suslov equations are obtained by re¬ 
duction from the nonholonomic equations produced by the Lagrange-d’Alembert 
principle. This particular reduction process leads to fewer degrees of freedom and 
also turns the dynamical equations from second-order differential algebraic to first- 
order. Moreover, the configuration manifold is a Lie algebra, concretely so(3), a fact 
which allows the direct application of many results on dynamical systems defined 
on linear spaces. Particularly, we focus on their discretizations, both direct ones 
and those coming from discretizations of the Lagrange-d’Alembert principle, and 
meticulously study the conditions under which the constraints are preserved. In 
terms of consistency of integrators, we wonder about the relationship between these 
in unreduced and reduced settings. We find out that this relationship is nontrivial, 
and can be completely worked out when we choose a particular retraction map re¬ 
lating 50(3) and so (3), say the Cayley map. The study of consistency also requires 
the introduction of a suitable metric on 50(3). The Suslov problem on 50(3), i.e. 
a rigid body with some of its body angular velocity components set equal to zero, 
is treated. General considerations and analysis of this system (and its generaliza¬ 
tions) may be found in EJ[I11E1]- Although we will not consider integrability issues 
(in that regard we refer to [lUEl]), we shall carefully study the consistency order 
of general and variational discretizations, and show their performance by several 
simulations (for extensive simulations by other reduced variational integrators we 
refer to m)- Finally, we address the discrepancies between continuous and dis¬ 
crete dynamics from a non-autonomous perspective, thereby taking into account 
previous results in Eziiin]. 


The paper is structured as follows: In ^ and we provide a comprehen¬ 
sive introduction to nonholonomic mechanics, variational nonholonomic integra¬ 
tors, symmetries in nonholonomic Lagrangian systems, the Euler-Poincare-Suslov 
equations and, finally, variational integrators for reduced systems. We prove the 
equivalence of the reaction forces, say Lagrange multipliers, coming from the unre¬ 
duced and reduced settings, see proposition |4.2[ Moreover, we provide a new ver¬ 
sion of the discrete Euler-Poincare-Suslov equations, defined on the Lie algebra 
through the incorporation of a retraction map, in corollary 4.4 establishes the 
link between the results in [23 and the reduced framework. First, 15.3 accounts 




ON THE DISCRETIZATION OF THE EULER-POINCARE-SUSLOV EQUATIONS IN SO(3) 3 


for the relationship between the consistency order of the unreduced and reduced 
setups, which is nontrivial and can be completely worked out when the discrete 
reconstruction equation is given by a Cayley retraction map. Then, in subsection 
§5.4[ we study general discretizations of the Euler-Poincare-Suslov problem which 
preserve the nonholonomic constraints, as well as sufficient conditions in the case 


of variational integrators guaranteeing such a preservation (proposition 5.2). The 
section accounts for an exhaustive study of the Suslov problem on SO{3). We 
give consistency results both for general and variational discretizations, and illus¬ 
trate their performance through some plots. Finally, 0 is devoted to the study of 
distribution-preserving discretizations of the Euler-Poincare-Suslov problem viewed 
as perturbation of the continuous dynamics. For this, we apply the result by Fiedler 
and Scheurle m to the system under consideration, and establish a bound of the 
perturbation produced in the unreduced dynamics by a reduced integrator, and 
conversely. 

Throughout the paper, we use Einstein’s convention for the summation over 
repeated indices unless the opposite is stated. 


2. Preliminaries 

2.1. The nonholonomic setting. We define a nonholonomic Lagrangian system 
as a triple {G,L,D), where G is a finite dimensional manifold, L : TG —>■ R 
is the Lagrangian function and D C TG is a constraint distribution^ which we 
assume to be a constant rank linear vector subbundle of TG and non-integrable 
in the Frobenious sense. In this work we are interested in the particular instance 
G = S'0(3); however, the results introduced in the preliminary sections refer to any 
Lie group, and therefore we stick to the notation G. Locally, the linear constraints 
are written as follows: 

( 5 , 5 ) = = 0 , l<a<m, (1) 

where i = l,...,n, are coordinates of TG and rank(T)) = n — m, with 

m < n. The annihilator D° is locally given by 

D° = span = ^“(g) dg'", 1 < a < m} , 

where the one-forms /i“ are linearly independent. 

In addition to the distribution, we need to specify the dynamical evolution of 
the system through the Lagrangian function. In nonholonomic mechanics, the 
procedure leading from the Newtonian point of view to the Lagrangian one is given 
by the Lagrange-d’Alembert principle. This principle says that a curve g : I G 
E, —Q describes an admissible motion of the system if 

^ [ L{g{t) ,g(t))dt = 0 

Jti 

with respect to all variations such that Sg{t) € Dg(t)i ti < t < t 2 and the fixed 
endpoint condition is satisfied, and if the velocity of the motion itself satisfies the 
constraints. It is to be noted that the Lagrange-d’Alembert principle is not varia¬ 
tional since we are imposing the constraints on the curve after the extremization. 
Thus, one may consider the intrinsic data defining the Lagrangian nonholonomic 
problem to be given by the triple (G, L, D). Using Lagrange-d’Alembert’s principle, 
we arrive at the nonholonomic equations, which in coordinates read 


d 

dt 


w) 


dL 


dg^ 

dTi9)g" = 0 , 


= Aa nfig), 


(2a) 

(2b) 
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where Aq, a = are “Lagrange multipliers”. The right-hand side of equa¬ 

tion (2a) represents the reaction forces due to the constraints, and equations (2b) 
represent the constraints themselves. For more details we refer to Hi. 

As mentioned above, the concern of this paper is the discretization of this kind 
of systems for G = 50(3) and, moreover, when they are left invariant under the 
same group. One of the advantages of this group, and all finite-dimensional ma¬ 
trix groups in general, is that they can be embedded in R", for n big enough. In 
[20] , the discretization of nonholonomic systems in R" was treated, more concretely 
discretizations that preserve the principal geometrical object, which is the distri¬ 
bution manifold D. Using previous results in isniiini, it was proven the existence 
and uniqueness of solutions for this problem, even in the non-autonomous case, 
by reducing the basic equations of motion ([^, say the nonholonomic DAE, to a 
set of ordinary differential equations, say the nonholonomic ODE, on D. Eurther- 
more, upon a main result in it is shown that any integrator that preserves 
the distribution may be understood as being generated by the exact evolution map 
of a time-periodic non-autonomous perturbation of the original continuous system. 
Thus, these results apply in the case of the considered matrix group, although some 
subtle issues appear when dealing with the group structure as it will be shown. In 
regards of the discretization, the following definition will be of main importance: 


Definition 2.1. By a {p, s) order discretization of the nonholonomic problem we 
understand a sequence of points {{pk^Vg^, Xk)} G k = 0,1,N — 1, N, s.t. 

Vgk G Dgk (in other words p°‘{gk)vgk = 0 for any k), and, moreover, \g{tk + e) — 
gk+i \ ~ \vg{tk + e) — ~ with min(r,/) = p and, moreover, 

\X{tk + e) — Afc+i| ~ 0(e®+^) with s > 0. 


By I • I we understand the usual Euclidean metric. We note that the notation 
\g(tk + e) — gk+i \ might be easily misunderstood, since in general G is not a linear 
space. However, the particular meaning of this metric on 50(3) will be clarified 
when considering the consistency of reduced systems. On the other hand, the 
dynamics of {g(t),Vg{t), X{t)) is determined by (i), after introducing the second 
order condition g = Vg and projecting the resulting differential algebraic equations 
onto D, which provides an ODE evolving on D and a unique in terms of the 
other variables. We note that the evolution prescribed by ([^ is local. 


2.2. Discrete mechanics and variational integrators. The variational inte¬ 
grators are a kind of geometric integrators for the Euler-Lagrange equations which 
retain their variational character and also, as a consequence, some of the main 
geometric properties of the continuous system, such as symplecticity and momen¬ 
tum conservation (see [laiMlISH]). In the following we will summarize the main 
features of this type of geometric integrators. A discrete Lagrangian is a map 
Ld : G X G —>■ R, which may be considered as an approximation of the action 
integral defined by a continuous Lagrangian L : TG —?► R, that is 

Ld{go,gi)^ L{g{t),g{t)) dt, (3) 

Jto 

where g{t) is a solution of the Euler-Lagrange equations, say 


d /dL\ dL 
dt \dg ) dg 


( 4 ) 


joining g{to) = go and g(to -|- e) = gi for small enough e > 0, where e will be 
considered as the step size of the integrator. 
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Define the action sum Sd ■ —>■ R corresponding to the Lagrangian Ld by 

N-l 

Sd='^ Ld{gk,9k+i), ( 5 ) 

k=0 

where G G ioi 0 < k < N, N is the number of discretization steps. The discrete 
variational principle states that the solutions of the discrete system determined by 
Ld must extremize the action sum given fixed endpoints go and g^- By extremizing 
Sd over gk, 1 < k < N — 1, we obtain the system of difference equations 

DiLdigk, gk+i) + D 2 Ld{gk-i,gk) = 0 , ( 6 ) 

or, in coordinates, 

+ ^{gk-i,gk) = 0 , 

where I < i < n, 1 < k < N — 1, and x and y represent the n first and n last 
variables, respectively, of the function Ld- 

These equations are usually called the discrete Euler-Lagrange equations. Under 
some regularity hypotheses (the matrix Di 2 Ld{gk^ fffe+i) is supposed to be regular), 
it is possible to define from a (local) discrete flow map : G x G ^ G x G, 
by FL^{gk-i,gk) = {gk,gk+i)- We will refer to the Fl^ flow, and also (with some 
abuse of notation) to the equations as a variational integrator. This kind 
of integrators are symplectic-momentum preserving in some sense, which makes 
them interesting from the Geometric Integration perspective (see [HISS] for more 
details.) 


2.3. Nonholonomic integrators. Discretizations of the Lagrange-d’Alembert 
principle for Lagrangian systems with nonholonomic constraints have been intro¬ 
duced in [ini |34j as a nonholonomic extension of variational integrators. To define 
a discrete nonholonomic system providing a discrete flow on a submanifold of G x G 
and, as well, a corresponding version of discrete Euler-Lagrange equations (§ , one 
needs three ingredients: a discrete Lagrangian, a constraint distribution D C TG 
and a discrete constraint space Dd C G x G. 


Definition 2.2. A discrete nonholonomic system is given by the triple (GxG, Ld, Dd) 

(1) Dd is a submanifold of G x G of dimension 2n — m with the additional 
property that 

Id = {{g,g)\g e G} c Dd. 

(2) Ld '. G X G ^ a is the discrete Lagrangian. 


We define the discrete Lagrange-d’Alembert principle (DLA) to be the extrem- 
ization of the action sum 

N-l 

^d= Ld{gk,9k+i) (7) 


among all sequences of points {gfelo-Ar with given fixed end points go,gN, where the 
variations must satisfy Sgk G Dg^ (in other words 6gk G ker pA) and {gk,9k+i) G Dd 
for all k G {0,..., N — 1}. This leads to the following set of discrete nonholonomic 
equations 


Dl Ld {gk , gk+l) T D2Ld(^gk—l, gk) — 9 {gk) , 

{gk,gk+i) G Dd. 


(8a) 

(8b) 


For the sake of clarity, the condition (8bI may be rewritten as 4>f{gk,gk+i) = 0, 
where (p’f : G x G —R is the set of m functions whose annihilation defines Dd 
(in other words, a suitable discretization of the nonholonomic constraints (§). 
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Equations ([^, where Aq = (A/c-|_i)q is chosen appropriately by projecting onto Dd, 
define a local discrete nonholonomic flow map Ff’f : Dd —t Dd given by 

Ph^idk-i^gk) = {9k,gk+i), ( 9 ) 

where gk+i satisfyes Q provided that {gk-i,gk) & Dd, if and only if a regularity 
condition, say the matrix 

DiD2Ld{gk,gk+i) Ai“(5/c) 

^2(/'d(5fc,5fc+i) 0 

is non-degenerate, is fullfiled for each {gk,gk+i) in a neighborhood of the diagonal 
of G X G. 

To obtain distribution preserving integrators within this setting is not straight¬ 
forward, since the equations Q are defined on the discrete space G x G and more¬ 
over determine a discrete flow on the discrete distribution Dd- Therefore, to locally 
relate {gk,gk+i) G G x G and {,gk,Vg,f) G TG is in order, a task accomplished in 
pn] by means of discrete difference maps p [S3]. 

Definition 2.3. A finite difference map p is a diffeomorphism p : U{Id) —t V{Z), 
where U{Id) is a neighborhood of the diagonal Id in G x G, and V(Z) denotes a 
neighborhood of the zero section of TG, i.e Z : G ^ TG s.t. Z(g) = Qg G TgG, 
which satisfies the following conditions: 

(1) p{Id) = Z, 

(2) To o p{U{Id)) = G, where tq : TG -G G is the canonical projection, 

(3) To o p\j^ = = 7r2|/^, where are the projections from G x G to G. 


Using them, the so-called velocity nonholonomic integrators are: 

h,:=poFffop-\ ( 10 ) 

defining a flow : Tg^^G -G- Ty^^^G; {gk,Vg^) ^ {h+i,Vg,+,)- In general, this 
flow is not D—presrving, but sufficient conditions for that to hold are given, which 
in general involve an appropriate redefinition of the discrete nodes {gk} '-G {gk} 
and the particular discrete constraints given by The result is stated in the 
following proposition (see [10]): 


Proposition 2.1. Assume that {gk,Vgfi) G If Dd is defined by 4>^ := pL°‘ o p : 
GxG —>■ IR, then (fffe+i, ) defined by the velocity nonholonomic integrator (101, 
i.e. Fkj^, belongs to . In other words, —?> 


3. Symmetries 

3.1. Symmetries of Lagrangian systems. We say that the Lagrangian is in¬ 
variant under a group action $ : G x G —)■ G, if L is invariant under the lifted 
action of G on TG, i.e. L o T$g = L. Such a symmetry allows to define a re¬ 
duced Lagrangian system on the reduced phase space TG/G = (G x g)/G = g, say 
F g —)■ E, (where we have employed the left trivialization mapping tr : TG ^ G x g, 
Vg (g, Tglg-i{vgf), where £ : G x G ^ G-, {g,h) ^ gh] Ig : G ^ G, h ^ gh, 
is the left action and g is the Lie algebra of G). More concretely, the reduced 
Lagrangian is obtained using the symmetry by 

L{g, g) = L{g-^ g, g-^ g) = T(e, ^ =: KO, (H) 

where we are employing the shorthand notation Tgikg '.= hg (we shall equally 
employ the notation Tgih = {Ih)*), e G G \s the identity element and g~^ g '■= f, G Q 
is called the reconstruction equation. Oviously, employing the right-invariance leads 
to equivalent results. 
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The reduction process allows to diminish the degrees of freedom of the given 
system and to obtain a first-order ODE describing the dynamics, instead of a 
second-order one. The paradigmatic example in the context of this work is the 
Euler-Poincare equations. The Euler-Lagrange equations Q are second-order, and 
obtained by the extremization of the action functional L{g{t), g{t))dt with fixed 
endpints Sg{ti) = 0, Sg{t 2 ) = 0 and arbitrary variations 6g G TgG. Now, we 
assume a left-symmetric Lagrangian. Thus, we consider the left-translated vari¬ 
ations Sg G TgG, namely g~^Sg = g, where g G Q. The endpoint condition 
6g{ti) = Sg{t 2 ) = 0 implies directly that g{ti) = 17 (^ 2 ) = 0. Furthermore, the 
variations 6^ are not free anymore but determined by the equations f = g~^ g and 
g~^ 6g = g. After an easy calculation we arrive at 


5f = g + ad^g, 


( 12 ) 


where the adjoint action on the algebra stands for ad^g = [^, 17 ]; [■, ■] : 0 x g —>• g 
represents the Lie algebra bracket. Now, the extremization of the action functional 

siO=rim)dt, (13) 

Jtl 

for an arbitrary curve g{t) in the Lie algebra with fixed endpoints g(ti) = 17 (^ 2 ) = 0; 
provides the so-called Euler-Poincare equations 


d /dl\ 
dt \ 9^ / 



(14) 


Under regularity conditions, these equations provide the solution curve f,(t) C g. 
Finally, the reconstruction equation g = g f will provide the configuration curve 
g{t). Thus, in case of a symmetric Lagrangian, the second-order equations Q may 
be replaced by the two first-order ones 


dt V 9^ / 



and 9 = gt 


which can be solved consecutively. However, they are not “Lagrangian” in the strict 
sense, since the variations (121 are not free anymore. 


3.2. The nonholonomic case: the Euler-Poincare-Suslov equations. In the 

nonholonomic case, besides a left-invariant L we shall consider as well a left- 
invariant distribution D (which we will call a left-left or 11 system). D C TG is 
left-invariant if and only if there exists a subspace g° C g such that Dg = {ig)*Q^ C 
TgG for any g G G. Let a“ G g*, a = 1, be independent annihilators of the 

subspace g°, i.e. 

g'" = {? G 0 I (a“, 0 = 0 : a = ■■■: W ■ (15) 

Consequently, the left-invariant constraints on TG are defined by the equations 
(a“ , g~^g) = 0. This last equation establishes the relationship between the non¬ 
holonomic constraints of a left-left system and the usual linear ones ([^. More 
concretely, the nonholonomic constraints are locally given by (j)°‘{g,g) = gT{g)g, 
and therefore 

(a “,0 = {a°‘,9~^g) = = (^*- 10 “, 5 ) = {9°‘{g),9) (16) 

Thus, fi°‘{g) = £*_!«“. The left invariance is easy to prove through 

<('“( 55 , 55 ) = (Ai“(55),55) = (^(3g)-ia“, (4)*5) 

= (^;_,£;_,a“,(4)*<?) = (^;_,a“,(V0*(4)*5) = {i;-.a^,g)=r{9,9)- 
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The dynamics of a 11 system is determined by the so-called Euler-Poincare-Suslov 
equations according to [25] : 


d fdl\ 
dt \ 9^ / 


= adj 


(a“,O=0. 



+ Ao 


(17) 


Furthermore, as in the unconstrained case, the dynamics of the group variables g is 
obtained through the reconstruction equation g = g^. We note that the equations 
in ( [T7| are the analogue to those in (H) in the case of 11 reduced systems. For more 
details on the Euler-Poincare-Suslov equations and their integrability we refer to 

mm- 


4. Variational integrators and reduced systems 


When one intends to discretize a Lagrangian problem defined on a Lie group, the 
most natural way is to directly discretize the equations Q . Through the discretiza¬ 
tion of variational principles, Marsden and collaborators obained in |28l [29] the 
discrete Euler-Poincare equations, which define an algorithm on the reduced space 
that is shown to be equivalent to the discrete Euler-Lagrange equations [531ISH] in 
the sense of reconstruction. This algorithm is obtained by assuming left-invariance 
of the discrete Lagrangian function Ld introduced in §2.2| More concretely, with 
respect to the induced group action on G x G we assume 


Ldi^gki gk+i) — Ld{g gkj g gk+i) 

for all 5 , gk, fffc+i e G. If we set 5 = g'^^, we get Ld{gk,gk+i) = Ld{gl^ gk,gl^ gk+i) 
= Ld{e,g^"^ 9k+i), leading to the quotient map 

TT : G X G ^ {G X G)/G = G, {gk,gk+i) 9k^ 9k+i- ( 18 ) 

This projection map defines the reduced discrete Lagrangian : G —>■ R for any 
G—invariant Ld hy Id o n = Ld', furthermore if we define the points in the quotient 
space as Wk ■= g^^ gk-eit we obtain 

^tl(bF/c) — ^didk 9k+l^ ~ ^dip, g^^ 9k+l) ~ ^d{gk, gk+i) ■, (1^) 


and the reduced action sum 

N-l 

Sd=Y. ^diWk). ( 20 ) 

k=0 

Applying calculus of variations to this action functional we obtain the so-called 
discrete Euler-Poincare equations ([31]), for the sequences {Wk}Q.]y_i, 

such that Wk = gk^gk+i- 

- + r^/diWk) = 0, k = 0...,N-2, (21) 


where iw^ and rw^+i are the left- and right-translation of the Lie group and ' 
denotes the derivative. 

Roughly speaking, the discrete Euler-Poincare equations (21 1 generate, under 
regularity conditions, a discrete flow Wk >->■ Wk+i on G, generating the sequence 
{lAfc}o-A-i> which, after reconstruction gk+i = gkWk, approximates the original 
flow g{t) defined by the original Euler-Lagrange equations (|^. To see the relation¬ 
ship between the continuous and discrete Euler-Poincare equations, (141 and (21) 
respectively, is not straightforward (one might argue, for instance, that they are 
defined in different spaces, g and G respectively) and needs more elaboration as 
shown next. 
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4.1. Discretization using natural charts. The above mentioned relationship 
between G and g is achieved by means of the so-called retraction map t : g — >■ G, 
which is an analytic local diffeomorphism near the identity such that r(^)T(—= 
e, where ^ S g. In this sense, r provides a natural chart on G, and Wk are 
regarded as small displacements on the Lie group, linking gk and gu+i- Thus, it is 
possible to express each term through a Lie algebra element that can be regarded 
as the averaged velocity of this displacement. Two standard choices for r are the 
exponential map and the Cayley map. Regarding ^ as the “velocity” of a given 
displacement in the group, it can be expressed as 

= 'r-^{9k^9k+i)/e = T-^{Wk)le. 

The difference Wk = 9 ^^ 9k+i G G, which in general is an element of a nonlinear 
space, can now be represented by the vector ^k- Moreover, the discrete Euler- 
Poincare equations can be expressed in terms of the velocities the retraction 
map and its derivatives [7], which are defined by: 


Definition 4.1. Given a retraction map t : q ^ G, its right trivialized tangent 
dr^ : g —> g cind its inverse dr^^ : 0 0; o,re such that for g = t(^) S G and 

r],^ G 0, the following holds 


= dr.t77r(^), 

diT-\g)g = dT^^{gT{-^)). 


Using these definitions, variations <5^ and 6 g are constrained by 

dfk = dr^j^(—-|-AdT(£^j,)?7fc-|-i)/e, (22) 


where g ^ rjk '■= 9k^^9k', this expression is obtained by straightforward differen¬ 
tiation of ^k = gk+i)/e. Note that gk forms the sequence {?7fe}o-Ar where 

Vo = Vn = 0 since Sgo = Sg^ = 0. Now, let us define a Lagrangian function 
Id ■ 3 ^ ^ SiS a suitable approximation of the reduced action functional (13). 


Remark 4.2. The naming of : 0 —>■ R is controversial. We are allocating the 
name discrete for : G x G —)■ R and : G —>■ R, the latter in the reduced case. 
These are considered as the discrete counterparts of the continuous Lagrangians 
L : TG —R and Z : g —R. We note that TG and g are the Lie algebroids 
associated to the Lie groupoids G x G and G, respectively. Thus, in the general 
framework introduced in [JS], we associate the discrete Lagrangian Ld : S —R to 
the continuous one £ : AS —>■ R, where 9 and AS are the linked Lie groupoid and 
Lie algebroid, respectively. Therefore, to call Zd : g —>■ R a “discrete” Lagrangian 
function could be misleading. 


After the introduction of these ingredients, the following result establishes an 
alternate version of the Euler Poincare equations on the Lie algebra. 

Proposition 4.1. Let Z^ : g —R 6e a Lagrangian function, a sequence 

in g s.t. fk = '''~^{9k^9k+i)/^! with a given retraction map t, and a sequence 
where gk = 9k^^9k- Then, the following assertions are equivalent: 

(1) {Cfelo N-i extremizes the discrete action 

N-l 

^ (23) 


for constrained variations Sf,k (22) and an arbitrary sequence {gk}Q.N 
go = r]N = 0. 
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(2) The discrete Euler-Poincare equations defined on the Lie algebra hold 
Ua) - mk+i) = 0, fc = 0, N-2. 

The last equation may be rewritten as 

(drr^)* f'(a) - m^+i) = 0, (24) 

where it has been taken into account that (dT^j^_J* = (dT“g^^_^)* (see 

m for the proof). This equations, which can be undersood as a variational in¬ 
tegrator for a reduced Lagrangian system, define, under regularity conditions, a 
discrete flow on the Lie algebra t Cfc+i; generating a sequence {Cfc}o- 7 v- 2 > 'which 
approximates the continuous solution of the Euler-Poincare equations (14). The 
discrete configuration variables are recovered by the discrete reconstruction equa¬ 
tion, namely 

gk+i = gkWk = gkT{ef,k). (25) 


4.2. The nonholonomic setting. The extension of the previous results, and also 
of those in [5], to the nonholonomic left-left setup defined on Lie groups has been 
developed in m- First, if the continuous distribution D C TG is left-invariant, 
it is natural to require that the discrete constraint submanifold C G x G is 
left-invariant with respect to the diagonal action of G on G x G, namely 

(l>d{9k,gk+i) = (I)di39k,ggk+i), Vg G G, a = 1, ...,m. (26) 

Setting g = g^^ , one arrives at the conclusion that there exist functions (pf ■ G ^ 
IR, q; = 1, ...,m such that 

<l^d{9k^ 9k,gk^ gk+i) = ipd{9k^9k+i) = <Pd(bPfe), (27) 

whose annihilation defines completely Dd- The main result is stated in the next 
theorem [15] . 


Theorem 4.3. Let : G x G —>■ R 6e a left-invariant discrete Lagrangian, Id : 
G — >■ R 6e the reduced discrete Lagrangian, and D C TG, Dd G G x G be the 
constraint distribution (also left-invariant) and the discrete constraint submanifold 
respectively. Then the following assertions are equivalent: 

(1) {gfcjg.jv critical point of the action ([^ for constrained variations Sgk G 

s.t. Sgo = SgN = 0. 

(2) {gfclg.jv satisfies the discrete nonholonomic equations 

D 2 Ld{.9k—It 9k) 4“ DiLd{gkt 9k+i) — 3 {9k)t (28a) 

<)d{9k,9k+i) =0, (28b) 


(3) 

(4) 


for k = 1, ...,N — 1 ,where Xa = (Afe+i)^ are chosen appropriately (this is, 
by projecting onto Dd-) 

{Wfclo-jv-i *■5 ® critical point of the reduced action sum (20), with respect 
to variations 5Wk induced by the constraint variations Sg^ G Dg,.. 
{Wfc}g.jY_i satisfies the reduced nonholonomic equations or discrete Euler- 
Poincare-Suslov equations; 


~''^Wk+Jd{^k+l) + iwjd{^k) — Xa a“, 

P^{Wk+i) = 0 , 


(29a) 

(29b) 


for fc = 0,..., N — 2, where again Xa = (Afc+i)Q are chosen appropriately. 


The constraints : G —>■ R determine the so called constrained displacement 
subvariety s na, say 


§ = {tTGG| p‘f{W) =Q) gG, 


( 30 ) 
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which completely determines Dd by Dd = {{g,g) £ G x G \ g G §}. 

With regard to the Lagrange multipliers in equations (281 and (p^ , to determine 


if they are equal needs more ellaboration. The right hand sides of ( |28a[ ) and (29a) 
represent the reaction forces provoked by the nonholonomic constraints, forces that 
involve the Lagrange multipliers. Thus, it is interesting to study whether in the 
unreduced and reduced setups these forces are equal. For that we need to check 
both if, in the two setups, the discrete dynamical equations and the constraint 
subspaces are equivalent. The first item is treated in the next result. 


Proposition 4.2. Equations of (28a) and (29a) are equivalent. 


Proof. We prove this in the unreduced-to-reduced way and by direct computations. 
By the left invariance of Ld we have 

7V-1 N-l 

X] Ld{gk,gk+i) = X! Ld{ggk,ggk+i)- (31) 

/c=0 


Taking variations of the right hand side we obtain 

N-l 

^ {DiLd{ggk,ggk+i), S{ggk)) + {D 2 Ld{ggk,ggk+i), 5{ggk+i )); 

/c^O 

furthermore, setting g = it follows that S{g'^^gk) = 5e = 0 and 


N-l 

{D2Ld{e,g^^ gk^i), Sg'^^ gk-^i + 5 ^^ Sgt-^i) 

k^O 

N-l 

= -9k^ 5gk g^^ gk-^i + 9k^ 

N-l 

= X! {^d(Wk), -{rwj*{£g-^)*5gk + {(-Wk)* , 

fc=0 


where in the last line we have taken into account that Wk = g^'^gk+i and ld{Wk) = 
Ld{e, Wk). Rearranging the summation index after setting 6go = SgN = 0 we arrive 
at 

N-l 

(^;-i i-r*wMWk)+iw/d{Wk-i)) , 6gk) , 

A.-1 

which, considering that Sgk € Dg ^, yields 


{-r*^J'd{Wk) + twJ'AWk-i)) = 

and consequently ( 29a[ ) just by shifting k i-A k + \. Taking variations of the left 
hand side of (31) we arrive directly at (28a), proving the claim. □ 


On the other hand, we note that the subspaces determined by (28b) and (29b), 
say Dd C G X G and § C G respectively, despite related to each other, are different 
in principle. Consequently we conclude that the multipliers resulting from the 
projection of the equations (28a) and (29a I onto these spaces will be different in 
general. 

Considering a retraction map and equation (24), the following result is a corollary 
of theorem (4.3). 


Corollary 4.4. Let 0 : 0 — t R and : g —>■ R &e defined as above and r a 
retraction map. Then, the following statements are equivalent: 
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(1) {Cfelo-AT-i ® critical point of the action functional (23) for constrained 
variations Sf,k (22) and a restricted sequence {rik}Q.j^ G 0^ s.t. rjo = tjn = 

0 . 

( 2 ) {Cfelo .pf_i satisfy the reduced nonholonomic equations or discrete Euler- 
Poincare-Suslov equations defined on the Lie algebra: 

^d(a+i) = o, 

for k = 0, ...,7V — 2 and where Aq, = (Afc+i)^ are chosen appropriately. 


Proof. By direct computations we obtain 

N-l N-1 

<5 H , dr-i(-?7fe + Adr(,^,)Vk+i)/e^ 

k^O k^O 

N-l 

= U^k-i) - {dr-^lr mk ), m/e) = 0 , 

k=l 

where in the last line we have rearranged the summation index taking into ac¬ 
count that m = Vn = 0. Using again that (dT^’^^_ J* = (drUjj^_^)* and 

considering that m ^ ^ arrive at 

(dr-i^J* Uik) - mk+i) = Aaa“, k = 0,..., TV - 2, 

where we have used the shift k ^ k + and therefore the claim holds. □ 


Note that the second equation in ( [3^ defines a subset, which we will denote 
by C 0 , 0 ^' := ker(^J^, such that in general 0 ^' 0 ^. Note as well, that the 

suitable Lagrange multipliers (Afc+i)^ must be chosen appropriately in this case by 
projecting onto 0 °'. The equations ( [^ define a discrete flow : 0 ^= —)• 0 ^=, 
^fe+i = F^^{f^k), only under some regularity conditions, conditions which may be 
obtained using the implicit function theorem and locally ammounts to the invert- 
ibility of the following matrix 


+ (dr-iNmk+i) 



(33) 


Using coordinates we can write the upper-left term as 

did 
dHd 








(f/c+l), 

£i i^k+l)- 


Here, we are defining locally the pull-back of the trivialized inverse tangent retrac¬ 
tion map as a linear operator by := we shall see that these 

definitions are useful in the case of matrix groups. 


5. Application to S0{3) 

The group iS'0(3) is formed by 3 x 3 real orthogonal matrices such that their 
determinant is equal to one (with some abuse of naming, we will sometimes call 
these matrices orthonormal). It can be also interpreted as the -|-1-determinant 
connected component of 0(3), which is the group of rigid transformation in H^. 
From the physical point of view this is relevant since, among other reasons, it is 
the configuration manifold of the the Lagrangian and Hamiltonian formulations of 
the rigid body problem. 
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The corresponding Lie algebra so(3), which is the algebra of real 3x3 skew- 
symmetric matrices, is isomorphic to the Euclidean space through the isomor¬ 
phism • : —>• so(3), w !-)■ w, given by 

( 0 -UJ3 UJ2 \ 

013 0 —uji j G so(3) (34) 

—UJ2 ^ / 

for uj = {u!i,u! 2 ,uj 3 ) G R^. In this representation, the antisymmetric bracket oper¬ 
ation is the standard vector product in R^ (namely [?), /t] = 77 x /r for 77, G R^). 

Furthermore, in following sections the definition of a distance between group and 
algebra elements will be relevant, respectively. In the latter case, it is common to 
pick the Killing form, since it is invariant under all the automorphisms of so (3). 
Namely, 

= -^trace (^ 77 ) = ^trace^^T)), (35) 

where ^,77 G so(3) and ^77 represents the usual matrix product. As it is well- 
known, this bilinear form corresponds to the usual Euclidean product on R^, say 
for ^,77 G R^, we have = ^ ■ V = Thus, the distance 

between ^,77 G so(3) may be defined by 

k - »7l = \i -fl\ = ii -vd- (36) 


Concerning the Lie group, to pick a distance is a subtle issue, even more if one wants 
to stick close to Euclidean metrics. In the reduced setting, the update of the Lie 
group elements by the developed variational integrators is given by flfc+i = Ofc T(ec5) 
(25), where 17^ G SO{3) and uj G so(3). In case of T=exp, it is clear that flfc+i 
would also belong to SO{2>). Nevertheless, for later applications, we will pick a 
different retraction map, for instance a suitable truncation of the exponential. This 
choice in general prevents flfc+i to lie in SO{?>). Following this idea, we establish 
the self-distance in SO{3) as 


dist(n, n) = |j — 

where we use the Euclidean metric of matrices, induced by |np = irij 


(37) 
p. As it 

is natural, if S7 G SO{3), then dist(n, 17) = 0. In other words, what (37) measures 
is how far is a matrix from being orthonormal in terms of an Euclidean metric in 
a matrix space. This generalizes to 


dist(f7i, 172 ) = \I — 17 ^ 172 ! 


(38) 


measuring how far is f7i from being S72 within SO(3). It is easy to prove that 
dist: SO{3) x SO{3) —?► R is indeed a metric, and moreover allows to establish the 
following lemma, which is very useful for our purposes 

Lemma 5.1. Consider f7i,I72 G SO{3) such that dist(I7i, 172) = Then 


f7i - 172 = 

where H is a 3x3 real matrix defined by 9 real entries Hij = 0{e'''^^)hij{e), where 
at least one of the nine hij is independent of e. 


^In case of SO{3) matrices it might be more natural to consider the Frobenius matrix norm 
||f2l —r22||f' = 2tr(/—Q)(^r22). However, this metric only allows to establish bounds for the diagonal 
entries of the difference matrix I — Qf Q2 , which is not enough to establish the consistency bounds 
in the context of this work. 
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Proof. Let us define = H, then it follows directly that ^1 — ^2 = iliH. On 

the other hand, the condition dist(r2i,02) = 0(e’'+^) implies that \H\ = 0(e’'+^), 
and therefore 



Thus, all the matrix entries may be factorized as Hij = where at 

least one of the entries hij must be independent of e, since otherwise the square 
root would be of order 0(e’'+^) as e —>■ 0. □ 


5.1. The Cayley map in SO{3). Concerning the retraction map r, we employ 
the Cayley map since it is simple and also computational efficient. Say, it is suitable 
for iterative integration and optimization problems since its derivatives do not have 
any singularities that might otherwise cause difficulties for gradient-based methods 
[IS] . The Cayley map cay : g —)■ G is defined by cay(^) = (e — |)“^(e -I- |) and 
is valid for a general class of quadratic groups. The quadratic Lie groups are those 
defined as 

G = {r G GL{n, R) I Y^PY = P} , 

where P G GL(n,]R) is a given matrix (here, GL(n, R) denotes the general linear 
group of degree n; 0{n) or SO{n) are examples of quadratic Lie groups). The 
corresponding Lie algebra is 

g = {n G gZ(n, R) \ pn + n^p = 0 }. 


The right trivialized derivative and inverse of the Cayley map are defined by 

dcay^ ?7 = (e - r]{e+ 

dcay"^? 7 = (e - |), 


for ^,7? G g. Using the identification (34), its particular form in the case of SO{3) 
is 


cay(d)) = I + 


1 


1 + lfP 


LJ - 


CJ 

T 


(39) 


One of the most interesting properties of the Cayley map is that its image belongs 
50(3), i.e. cay(w) cay(d))^ = cay(w)'^ cay(a)) = I. Furthermore, the linear maps 
dT,t and dr^^ are represented by the 3x3 matrices 


dcay,^ = dcay^^ = I - Y + 


We point out to the reader that the same notation | • | will be used for the Euclidean 
norms on R^, so(3) and matrices. Taking into account (37), Ofc+i = Ofccay(a)) and 
(39), after a straighforward calculation we arrive at 

dist(Ofc+i,Ofc+i) = 0 (40) 


for Ofc G 50(3). Therefore, the update itk+i defined by (25) in the case r =cay 
remains in 50(3). 


5.2. The Suslov problem on 50(3). The Suslov problem is a well-known exam¬ 
ple of a left-left system, introduced in [41]. It describes the motion of a rigid body 
suspended at its center of mass in the presence of a constraint that forces the pro¬ 
jection of the body angular velocity along a direction fixed in the body to vanish. 
The generalization to the group SO(ri), called the generalized Suslov problem, has 
been described, both from continuous and discrete points of view, in danainiiiii. 
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Let I = (%) be the inertia tensor I : so(3) —>■ so(3)* of the body, 1“^ = 
(Jb) its inverse, and w G be the body angular velocity vector. The dynamics 
is determined through the reduced Lagrangian ( [IT| ), which in this case reads I : 
so(3) —t R 

= ^(Iw,w). (41) 

Let a be the direction, fixed in the body, of the vanishing component of the angular 
velocity. Thus, the nonholonomic constraints read (a, uj) = 0, where we are using 
the bracket as the pairing between (R^)* and R^. Taking into account ( [T7| ) and 
(41), the Euler-Poincare-Suslov equations, written in R^, are 


Iw = Ia;xa; + Aa, 

(o, uj) = 0, 

where A G R is the Lagrange multiplier. Componentwise, we have 

\ / Wa (l2iWi) — 012 (IsiWi) \ / «! 


I 


(42) 


UJ2 I ~ I *^1 ~ ^3 (IliWi) I + I 

<^3 / \ <j-'2 (IliUJi) — UJi / \ 03 

and oi wi + 02 012 + 03 W 3 = 0. Without loss of generality, we may choose a as the 
third component of the body frame { 61 , 62 , 63 } in R^, say a = (0,0,1). Then, the 
constraint becomes 013 = 0 and 

0 0 012 

so(3)^={| 0 0 -wi |}cso(3). 

—012 Wi 0 


Moreover, (42) becomes 

— UJ2 (IsiWi) 

Wl (S-SiUJi) 

UI2 (IliWi) — UJi (l2i(^i) 

where, from now on, i = {1,2}. After a straightforward computation we arrive at 




(43) 


Wl = — TT-T (I22W2 

I I 


1 


U}2 — 


|I,. 


Il2<j-’l) (IsiWi) , 
(I 21 CLI 2 + IllWi) (IsiWi) , 


(44) 


where Im = 


Ill 

I21 


I12 

I22 


A(w) = wi (12*011) - 02(111011 


is considered to be non-degenerate; moreover 
IszO'l 


|I. 


((I32I21 — l3ll22)o’2 + (I32I1I — l 3 lIl 2 )o'l) . 

(45) 


In other words, what happens to the Euler-Poincare-Suslov equations (42) in the 


presence of the nonholonomic constraint 013 = 0 is that they decouple into a dif¬ 


ferential part (44), this is a system of nonlinear ODEs which we will denote by 


o = /(o') for simplicity, and an algebraic part (45). 


5.3. Consistency order. Now, we extend the consistency results from [20] to the 
present case, when G = SO{3). We consider a general Euler-Poincare-Suslov prob¬ 
lem in SO{3) and therefore we may have two scalar constraints. Thus, we keep the 
a“ notation instead of the single a as in the previous subsection, where we intro¬ 
duced the specific example of the Suslov problem. In the definition |2.1| we have in¬ 
troduced the consistency order of a general integrator of the nonholonomic DAE p). 
However, the reduction process leads to the Euler-Poincare-Suslov equations (17), 


which are first-order instead of second-order equations. We take into account all 
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these facts to establish a relationship between the {p, s) consistency of a discretiza¬ 
tion of the nonholonomic DAE ([^ (for instance velocity nonholonomic integrators 
( [l0| )) and a general discretization of ( [T7| , and conversely (always in the left-left 
framework). For that purpose, we employ the notions of distance ( |36| ) and (37). 

implies |D(tfc-|-e) - flfe+il = |/- D(tfc + e)'^nfe+i| ~ 0(e’'^), 
~ 0(e*+^) and |A(tfc + e) — Afc+i| ~ 0(e®+^), where the 


Namely, definition 
\vn{tk + e) - 


2.1 


fc "T cj — Ak^l\ 

continuous dynamics (f2(t), ua(t), A(t)) is determined by ([^ and the discrete one 
(rife+i, , Afe+i) by the prescribed integrator. On the other hand, the notion of 
(p, s) order integrator may redefined for the Euler-Poincare-Suslov equations to be 
an integrator such that \ui{tk+e)—u)k+i\ ~ and |A(<fc-|-e) —A^+i| ~ 0(e®+^), 

where (a)(t),A(t)) is given by 0 and (wfe+i,Afe+i) by the prescribed integrator. 
The relationship between both the continuous dynamics is established by the re¬ 
construction equation, i.e. vn{t) = 0(t) = 0(t)<i(t), while the discrete ones are 
linked through Ofc+i = Of, cay(ea}/c). Furthermore, we note that the left translation 
of a) S so(3) by a group element O € SO{3), i.e. Ow, is nothing but the matrix 
product rio;, where G SO{3) and w G IR^ as established by the hat isomorphism 

The link between consistency orders at the group and algebra level, which is not 
straightforward, is discussed in the following result. 


Theorem 5.2. Consider a (p, s) integrtor of the nonholonomic DAE (§ and a 
{p,s) integrator of the Euler-Poincare-Suslov equations ( [T7| ), when G = 50(3), 
1{Cj) = |(Iw,a;), for ui G so(3), and such that flfe+i = f2fecay(ewfc). Then the 
following assertions hold: 

1. the {p,s) integrator (DAE) generates a (p, s) integrator of 0. 

2. the {p,s) integrator (EPS) generates at least a (min(p, l),s) integrator of 

<§• 

Proof. In first place we take into account the dynamical variables, say n(t), vaft) 
and ujft). 

1) By construction, the (p, s) integrator (definitio n |2.1[ ) implies \I — Dk-\-i^{tk + 
e)^| ^ and therefore it follows from lemma [5.1| that 

+ e) + D,(tk -\- e), (46) 

where we identify H directly with 0{e'^'^^). On the other hand, regarding the 
velocities we have VQ,{tk + e) — = 0(e*+^). (This last expression could be 

misleading from the geometric point of view, since we are adding vectors belonging 
to two different vector spaces, say Tf2p,,^£)50(3) and Tn,,^j50(3); however, from 
the reconstruction equation these are nothing but rotated vectors in E,^, i.e. D.{tk -\- 
e) uj(tk -h e) and Ofc+i cvk+i and therefore we can understand the sum as embedded 
in E^.) Thus, taking into account the reconstruction equation we arrive at + 
e) ~ 'I’Ofc+i = + e)u!(tk -he) — flk+icok-i-i = 0(e^~^^), where the right hand side 

is identified with a vector in E^. Furthermore 


w(tfe + c) ~ ^(tk + e)^^fc+i<j-’fc-i-i — -|- e)'^0(e‘'^^). 


Now, considering (46), fl(tk -h e) = fl(tk) -h e(l{tk) -h 0{e^), oj{tk -I- e) = uj{tk) -h 
eui(tk) -h O(e^) and fl{tk) = uj{tk) = Wfc, from the last expression we obtain 

w(tfe -l- e) — u)k+i = f4(t/j -|- e)^0(e’’^^) -|- D,{tk -h e)"^D,{tk -h e) uJk+i 

= D^O(e"+i) + nl 0(e'+i) uJk + 0{e^+^) + 0(e’'+2). 

Now, taking the norm on both sides and noting that p = min(r, 1), we finally obtain 

\Co{tk + e) - Wfc+i| = \uj{tk + e) - u)k+i\ ^ 0(6^+^). 
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2) First we consider the consistency order in the group, namely |/ — ^k+i^{tk + 
e)^|. For this, it will be enough to set up the Taylor expansion n{tk + e) = ^{tk) + 
eCl{tk) + to second order. On the other hand, Ofe+i = Ofccay(ecbfc) = 

flk (j + ea^LJk + Y where = 1/(1 + e^^) = 1 - +0(e^). Taking 

this into account and considering that = —uj for w G so (3) we have that 

nl^,=nl-eu;knl+^^u;lill + 0{e^). 

Moreover, since Ct = Ow, we have Ct{tk) = ft{tk) u}‘^{tk) + Ofcw(tfe). With all these 
ingredients, recalling that fl{tk) = 0^ and after some computations we arrive at 

£2 ^ 

+ e) = / + —ojitk) + O(e^). 

In general, uj(tk) ^0 and consequently, 

|/-0^+i0(4 + e)|~0(e2). (47) 

However, if uj{tk) = 0 we get \I — ilk+if^{tk +e)'^\ ~ O(e^) and then the consistency 
order is at least min(p, 2). Thus, we can say that the minimum consistency order 
is given by min(p, 1), and higher orders depend on nontrivial cancelations. 

Regarding the velocity, since we are considering a p—th order integrator in the 
algebra, i.e. u}{tk + e) — uik+i = 0(6^+^), the following chain of equalities holds 

oj(tk + e) — uJk+i = ^{tk + + c) oj{tk + e) — VLl.j^iVLk+iUJk+i 

= n{tk + tfvnitk + e) - = 0(6^+^). 


Thus, Hfc+iH(4+e)'^un(4+e)-Pnfe+i = ^k+i 0(eP+^) and furthermore, employing 
(47) and that ftk+i = cay(ewfe), we obtain 

vn{tk + e) - VQk+i = ^k cay(ewfc) 0(6^+^) - 0{e'^)vn{tk + e) 

= HfcO(eP+i) + 0(eP+2) - 0{e^)va{tk) + 0{e^). 


Taking norms, we arrive at \vnitk + e) — "cnk+i | ~ (The same 

argument presented above regarding the sum of vectors at TQ((j,_|_g)S'0(3) and 
Tq^^^SO{3) holds.) 


Finally, we consider the algebraic variable, i.e. the Lagrange multiplier A. Since 
we consider the left-left framework for SO{3), equation (16) reads 

(a“,w) = (a“,H^H) = (a“, (^a-)* 


which in other words means that the constraint distributions in both the unreduced 
and reduced setups are related by the group lift, determining the same subspace of 
so(3)° C so (3), a subset prescribing a particular subspace so* (3)° C so* (3) through 
the Killing metric. Thus, the reaction forces involving the continuous dynamics of 
X(t) will also be equal if the continuous dynamical equations are the same. This 
can be proved in the case l{ui) = \OLuj,uj) noting that the left hand side of the 
reduced and unreduced equations, both elements of so* (3), say 



_d /_^A 
dt \duj) 


c^\ _ 

My ~ My 



Aaa“, 

Aaa“, 


are equal, a fact that follows from straightforward computation setting L{^l, Ct) = 
Consequently, the consistency bound |A(tfc + e) — A/c+i| ^ 0(e'^+^) 
remains the same. This assures the complete claim. □ 




18 


F. JIMENEZ AND J. SCHEURLE 


Remark 5.3. We observe that the dynamical part of theorem 5.2 
p i-A p, p i-A min(p, 1 ) consistency relationships 


namely the 
is general for any kind of La- 
grangian functions, since in any case the unreduced and reduced continuous setups 
are related to each other by the reconstruction equation. In the case of the Lagrange 
multipliers, the consistency order depends on the particular Lagrangian function, 
and therefore the generality of the result is left as subject of further study. 


Remark 5.4. We note that the result in theorem |5.2| is generalizable to any group 
which can be embedded in R", for n big enough. However, when dealing with a 
general Lie group G the situation is more subtle, since to define a suitable distance is 
not easy in principle. In terms of Riemannian geometry, nevertheless, the distance 
is always prescribed once we have defined a positive definite bilinear form in the 
algebra (•,•)£ : fl x g —>• R+. Furthermore, there is a bijective correspondence 
between left-invariant metrics on a Lie group G and the inner products on the Lie 
algebra g (see [Ml ETj for further details). This leads to the definition of a left- 
invariant Riemannian metric on G (•, ■)g : TgGxTgG —> R_|. through left translation, 
say 

{u,v)g = ((4-i)*m, (4-i)*'u)e, 

for any u^v G TgG^ and furthermore to the Riemannian distance 


ft2 

b 


dist{g,g)= I ( 7 (t), 7 (t))^,^, dt. 


(48) 


where 7 : [^ 1 ,^ 2 ] —>■ G is the geodesic with j(ti) = g, 7 (t 2 ) = g, such that g,g are 
close enough in terms of natural charts (j ]4.1[). This fact allows to relate a (p, s) 
integrator for (|^ and a {q, s) integrator for (171 in both directions, say unreduced 
to reduced and conversely, by the following lemma. 


Lemma 5 . 5 . Let (•, Oe ■ 0 x g —?► R+ be a positive definite bilinear form, (•, ■)g : 
TgG X TgG —>■ R+ be the induced left-invariant Riemannian metric in G and 
dist(-, •) : G X G ^ R+ a Riemannian distanee on G as defined in (481. Given 
\f,{tk + e) — ^fc+i| = 0(e’'“''^) the geodesic distance between f,{tk e) and then 

dist(p(tfc -I- e),gk+i) = 0(e’’+^). 

Proof. Let g : [^ 1 ,^ 2 ] —>■ 0 be the geodesic curve such that g{ti) = f{tk + e) and 
vih) = ■ffc+i, dehned as 'y{t)~^j{t), where 7 : [^ 1 ,^ 2 ] —> G is the geodesic with 
7 (G) = g{tk + e), 7 (^ 2 ) = Sfc+i- Then 


pZ2 

dist{g{tk + e),gk+i) = / dt 

Jti 

= f iv{t),g{t))edt, 
Ju 


rt2 


(7(t) ^7(t),7(0 ^i{t))-,Mdt 


where we have used the left-invariance of (•, ■)g. Thus, the right hand side is by 
definition the geodesic distance between ^(tfc -|- e) and ^fe+i, leading to dist(p(tfe -I- 
e), Pfe+i) = 0 (e*+^) as claimed. □ 


However, this is not enough to establish the full relationship between the consis¬ 
tency of integrators, since the bound of \vg[tk + e) — I needs to be established 
also in both directions. This is left as subject of further study, but necessarily this 
bound has to depend on the structure of the bilinear form (•, -je and its link to the 
particular Riemannian metric (•, ■)g. 
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5.4. Preservation of the constraints under discretization. The left trivializa- 
tion is a powerful tool when understanding the preservation of the constraints under 
discretization in the reduced setting for Namely, let us consider the initial 

data (rikjtbk) satisfying the nonholonomic constraints (15), this is {a°‘,Cjk) = 0 
and thus ujk € so(3)^. Consequently, by left trivialization, the associated pair 
{Qk,^k) = (flfcjflfcWfc) also satisfies the original constraints, i.e. = 0. 

Thus, we can construct flfc+i through the retraction map as suggested by the re¬ 
duction setting, say flk+i = ^kT{tCjk)- Finally, we just need to define ujk+i in terms 
of the previous data such that it satisfies the constraints, i.e. CJk+i G so(3)^. This 
can be done by choosing a suitable discretization of namely: 

DREPS(wfc,(A„) fe+i; Wfc+i) — 0, 

(a“,Wfc+i)=0, 

where by DREPS(a}fc, (Aajfc+i; Wfc+i) = 0 we denote a discretization of the first 


equation in (17) (the initials are named after Discrete Reduced Euler-Poincare- 


Suslov equations). In general, this discretization of the Euler-Poincare-Suslov equa¬ 
tions can be understood as a mapping DREPS: so(3) x R™ x so(3) —)■ so(3)*. Eor 
instance, the variational procedure, through the corollary |4.4[ provides a particular 


DREPS(wfc, (Aa)fe+i;a}fc+i) = 0 by means of the first equation in (32), namely 


DREPS(c5fe, (A„)fe+i;d;fc+i) = (dr-i-J* I'A^k+i) - (A„)fc+ia“. 

(50) 


We generally denote the coupled discrete equations (49) by 


DREPS(wfe, (AQ,)fe+i;c5fc_|_i) — 0, 

and assume them to be regular enough to provide uJk+i G so(3)^ and (AQ,)fc+i 
in terms of uJk] particularly, this regularity condition may be described locally, 
according to the implicit function theorem, by the non-degeneracy of the matrix 


dj, DREPS(a}fc, (AQ,)fc+i; Wfe+i) di DREPS(a}fc, (AQ)fc+i; Wfc+i) 

a“ 0 


(51) 


for close enough and Wfc+i, where di denotes the partial derivative with re¬ 
spect to the Tth slot in DREPS(-, •; •). We note that (33) coincides with this 
last regularity condition when we choose (50) and ip'^i'^k+i) = (a“,a)fe+i), con¬ 
sequently V(p2(a)fc+i) = a“. This process defines a discrete local flow {Q.k,Cjk) t 
(Dfe+i, a)fc_|_i), and furthermore (D^, Dfc) i—)■ {flk+i, Dfc+i), which schematically leads 
to the following algorithm. 


Definition 5.6. 

(1) Input data {Q.k,0Jk) s.t. {a°‘,uik) = 0, 

(2) Set Itk — 

(3) Define Dfe+i = Dfc r(ea)fc), 

(4) Obtain ihk+i from DREPS(wfe, (AQ,)fe+i; Wfc+i) = 0, 

(5) Output data (Dfc_|_i,a}fc+i), 

( 6 ) Set Ctk+i = D.k+iilik+ 1 , Output data {flk+i,Clk+i)- 


Some remarks are in order. First, we note that the discrete flow (flk,uJk) >—>■ 
(Dfe+i,Wfc_|_i) is well-defined on SO{3) x so(3), also (Dfc,Dfc) i-A- (Dfe+i,Dfe+i) on 
TSO{3), for a general r and e small enough. However, since r is a diffeomorphism 
only at a local level, the algorithm |5.6| does not preserve the original distribution 
D C TSO{3) established by the unreduced constraints = 0, in spite 

DREPS(wfc, (AQ)fc+i; tD/c+i) = 0 does preserve the reduced constraints (a“, uJk+i) = 
0. In any case, as shown in (40) Dfc+i G SO{3) when r =cay, and thus these two 
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discrete flows remain in the group when we restric to the Cayley map. Regarding 
the preservation of the constraints, we prove the following result. 


Proposition 5.1. Consider the algorithm |5.d| such that t =cay and e is small 
enough. Then, the output pair (flfe+i, satisfies the unreduced constraints 

= 0 . 


Proof. We set = flkcay{eujk), equations (37) and (40) allow us to write 


^fc+i^fc+i ~ ^■ 

Now, considering the unreduced constraints and employing a matrix notation we 
obtain 

Thus, since we consider constraints-preserving integrators in the reduced equations 
we have that (a“, uik+i) = 0 and the claim holds. □ 

Second, it is interesting to note that the variational procedure described in the¬ 


orem and corollary 4.4 does not provide directly so(3)^—preserving integra¬ 


tors. More concretely, the discretization of the constraints provided by (32) reads 
(p‘f{uJk+i) = 0, generating a perturbed subset so(3)^ 


C SO(3) by kert^^. 


It fol¬ 


lows quite obviously that the so(3)preservation is obtained, and thus so(3)°' = 
so(3)°, if we choose : so(3) —>■ R as As an interesting feature, 

we also point out that there is no need in this case to redefine the nodes {wfc}o.Ar_i) 
as happens in general for nonholonomic systems defined in Q = R” (see [20] for 
further details). 

The preservation of the constraints in the variational setting, when we restric 


ourselves to the discrete Euler-Poincare-Suslov equations (theorem 4.3), is less ob¬ 
vious. Locally, these equations generate a flow inside the submanifold § C SO{3) 
(30) if a regularity condition is satisfied. Again, this condition is obtained employing 


the implicit function theorem and is equivalent to the nondegeneracy of the matrix 


"^Wk+i ff{Wk+i) + {Wrw,^,rUWk+i) 
yp^iWk+i) 


0 


(52) 


where [Vrwk+i )* makes sense locally, in analogy with the case of the Euler-Poincare- 
Suslov equations defined on the Lie algebra and their regularity condition (33). 
Moreover, the relationship between § and Dd C 5'0(3) x SO{i) is established 
through the reduction process, say Wk — fl^flk+i, where {Ttk,^k+i) G Dd C 
SO{3) X 5'0(3) if and only if Wk & S. Thus, the discrete Euler-Poincare-Suslov 
equations generate a local flow §—>■§, which at the same time defines a local flow 
within Dd. That this flow preserves the original constraints requires therefore a 
local relationship between SO{i) x S'0(3) and TSO{3), a relationship that we may 
setup, as introduced in j ]2.3[ through a finite difference map p (definition |2.3[ ). More 
concretely, p : SO{3) x SO(3) —>■ TSO{3) gives 

(Rfc 1 id/c+i)' t ). (53) 

Taking into account all these elements, we can establish at this point a result con¬ 


cerning the D—preservation of the discrete Euler-Poincare-Suslov equations (29). 


Proposition 5.2. Let : SO{3) x SO{3) —> R (8b) be defined by := o p, 
where p is given by (53) and are the original one-forms /i“(fl) = (fl^)*a“ 
considered as functions : TSO{3) —)■ R. Then, the discrete Euler-Poincare- 


Suslov equations (29) generate, under the regularity condition (52), a D—preserving 


sequence flfc+i). 














ON THE DISCRETIZATION OF THE EULER-POINCARE-SUSLOV EQUATIONS IN SO(3) 21 


Proof. The discrete Euler-Poincare-Suslov equations generate the sequence {VEfe}g.jY _2 
evolving on §, which by “unreduction” generate {(rifc, evolving on 


according to Wk = Thus, cj)°'{V,k,^k+i) — 0. Furthermore, (53) implies 

that 

(jp{p.k, f^fe+i) = = 0, 

and the claim follows. □ 


Note that (53) gives us some freedom concerning the choice of that is the 


nodes for discrete integration (a fact that also showed up in the case of R" treated 
in [10]). Nevertheless, the natural choice for p is 

It is easy to understand that this choice of a finite difference map establishes a 
local isomorphism § = so(3)° through the retraction map t, since T~^{iljilk+i) = 
T~^{Wk) G so(3)° for Wk G §. 

6. Discretizations of the Suslov Problem on 50(3) 

In this section we carefully study some particular discretizations of the Suslov 


problem in 50(3) introduced in (5.2 


6.1. General discretizations DREP^ It is interesting to note that any dis¬ 
cretization DREPS(a;*, Afc+i; w^''"^) = 
respecting the local regularity condition 


(j^ (50) (which implies = wf = 0), 

1 (|51|) and applied to (|43|, leads to partic¬ 


ular discretizations of (44) and (45). Conversely, a general discretization of (^44 


prescribes a DREPS of the Euler-Poincare-Suslov equation on 50(3), for which we 
can derive an interesting result. Prior to its statement, we recall that w = /(w) is a 
system of nonlinear ODE defined in a vector space, and consequently we can apply 
a numerical method of arbitrary consistency order (Euler, midpoint rule, Runge- 


Kutta, etc.). Furthermore, the decoupling of (43) into differential and algebraic 


parts allows some freedom in the choice of its DREPS. In particular, we choose 


^l(a;^a;'=+^e) 

^2(a;^a;'=+^e) 


(54) 


Afe+i — A(a;^~''^), 


where we pick Afc+i according to (45) as a natural choice. We assume li,l 2 to be 


smooth functions chosen in such a way that they generate a p—th order consistent 


one-step discretization of (44), i.e. \uj{tk + e) — ^ with p > 1. 

Namely 


, _ , ,k 1 

^^ = 7^ (I22^l(^c^cc^■+^e) -Il2^2(w^a;'=+^e)) , 


IR 


- Wn 


1 


(55) 


— - - — 77—I (—l2i^i(w*, e) -I- , e )) ■ 

^ l-^m I 

Now we have all the necessary ingredients to establish the result. 


Proposition 6.1. Any p—th order discretization (55), p > 1, of (44) generates a 
discretization DREPS(a;^, Afe+i;a;^+^) = 0 (|54|) of the Euler-Poincare-Suslov prob¬ 


lem defined on 50(3), of order (p,p) for (0, and consequently at least (min(p, l),p) 
for ([^. 

^In this subsection and the next one we are going to raise the index k in the uj variables, to 
avoid any misleading mixing with the index, say uii. 
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Proof. The dynamical part is obvious, namely the order \ ~ 

is given by assumption. On the other hand, given the choice A^+i = 
according to (451 and the uj consistency bound, using the Taylor expansion of A we 
have 


|Afe+i — A(tfc + e)| — |A(a;^“''^) — X{uj{tk + e))| 

= \X{uj{tk + e) + 0(e^^^)) — A(w(tfc + e))| 

= |0(e^’+i)VA(a;(4)) + 0{e^+^)\^ 0{e^+^). 


Once the {p,p) integrator is confirmed for the reduced system, the second statement, 
this is the (min(p, 1), 1) consistency for the unreduced problem, follows directly from 
Theorem 15.21 □ 


Remark 6.1. The previous result may be refined concerning the algebraic part, 
generating a (p, s) integrator with s > p. For this purpose, we set 

Afc+i = A(a;^'^^) + e) 


in (54), with l\ chosen as shown next. In the new scenario, we have (where we omit 
the l\ variables for the sake of simplicity) 


|Afc+i — X{tk + e)| — |A(a;*''’^) — X(uj{tk + e)) + ^a| 

= |A(w^+i) - A(a;'=+i + 0{eP+^)) + h\ 

= \0{tP+^)D A(a;'=+i) + 0{e^P+'^)D'^X{uj^+^) + 0{e^P+^)D^X{uj'^+^) + • • • + 0(e"+i) + /a|. 


where we take the Taylor expansion of A(a;^“''^ + 0{eP~^^)) up to the term 
(s = np + n — 1, with n an integer). Thus, it is obvious that choosing l\ such that 

C>(eP+i)VA(a;'=+i) + 0(e2^’+^)V2A(w'=+^) + 0(e3P+3)V3A(a;'=+i) + --- + 0(e®)+lA = 0, 


we obtain the (p, s) integrator as claimed. 


6.2. Variational nonholonomic integrators. As presented in f|^ the variational 
integrators setting provides us with a framework for the numerical integration of re¬ 
duced systems. In particular, corollary |4.4| prescribes a particular DREPS, namely, 


(56) 


Regarding the discrete constraints (pd, although we have some freedom for their 
choice, (pd{.‘) = {a, •) is the most suitable, since it implies the reduced constraints 
preservation. 

For the case concerning this work, we will assume furthermore that UI 3 = 0, 
which yields the following form of the right trivialized inverse derivative of the 
Cayley map in 50(3): 


^ + T^i ~2 

dcay7j = I ^ia;iW2 ^ 

|W2 


W2 


1 + 

-fwi 


2^1 

0 


(57) 


As pointed out in the discussion below the equation (33), we see that the right tri¬ 
vialized derivative of some retraction maps, in this case the Cayley map in 50(3) is 
given by a matrix, and therefore in local coordinates it has two indices: (dcay^^)^ ^. 
Setting Idipj) = eZ(cb) as a first order approximation of the action s(tD) = 
l{u;) dt (13), where I : so(3) —>■ R is given by (41), and applying ( |^ with 







ON THE DISCRETIZATION OF THE EULER-POINCARE-SUSLOV EQUATIONS IN SO(3) 23 


(pd{,‘) = (oj •) (note that the constraint is single in the Suslov problem), we obtain 
the following algorithm 


I 


m 


,fc+l 

,fc+l 




—UJ 


,/c+l 

-^2 

fc + 1/ 


(I 3 


, ,k+l 
fc+1 


) + <^2 (Isiwf ) 






to 


fc+l, ,/c+l 


Wn 


, ,fc+l fc+1 

CJ-| UJcy 

k+l\2 




(CJ: 


fc\2 
1 ) 

,k. .k 
■>1^2 


k, ,k 




1^2 

,k\2 

^ 2 ) 




UJ, 



(58) 


which accounts for the discretization of the dynamical part of the Euler-Poincare- 

Suslov equations (we recall that 1^=1 )), while 

V 1^21 1^22 / 


Afc+i — 2 


fc +1 

1 


(l2,a;+i) +ccf(l2.+)) - - (a;+i(lHCc+i) +..2^(Ii,a;0) , 


(59) 


where the rescaling A^+i i-+ — Afe+i/e^ has been introduced, accounts for the alge¬ 
braic part (note that in this case we do not have freedom to fix the algebraic part 
since it is directly given by the variational integrators setup). This rescaling may 
be understood in the context of the construction of the variational integrator from 
a discretization map p : SO{2,) x 5'0(3) —>■ TSO{3). More concretely, as can be 
shown in view of proposition |5.2[ any discretization (p°‘ (x p°‘ o p preserves the con¬ 
straint for the discrete flow (17*,, 17^) i-+ (17^+1, f7fc+i). Thus, setting <))“ = op 

accounts for the mentioned rescaling. 

Once we have derived the discrete equations, we may establish the following 
result regarding their consistency order with respect the continuous dynamics. 


Proposition 6.2. The numerical method (58), (59), is a (2,*) consistent method 
with respect to the Euler-Poincare-Suslov problem (17) defined on 50(3), and con¬ 
sequently (1,*) with respect to ([^; where by * we mean that the integrator is not 
neccessarily consistent with respect to the multipliers. 


Proof. To prove this result, we just consider the Taylor expansion uj(tk -|- e) = 
uj{tk) euj{tk) ^uj{tk) 0{e^), where ti; is given in (44), and compare, order by 
order, with provided by ( [^ . In first place, we calculate the second order 

time derivatives, namely 


wi = - 


1 


lIrnP 

1 

++ 

1 


{l22(l3iWi)^ (121^2 -b IllWl) — Il2 (Isi^i)^ (I 22 W 2 4- Il2+’l)} 

I 22 1+2 (Isi+i) {~l3l(l22+2 + Il2+l) + l32(l22+2 + Il2+l)} 

I12 +i(]l 3 i+i) {—l 3 l(l 22+2 + I12+1) + 132(122+2 + I12+1)} , 


and 


+2 


——^ {l2l(l3i+i)^ (121+2 + Ill+l) — Ill(l3i+i)^(l22+2 + Il2+l)} 

l-^m I 

+ 7 +—I 21 +2 (l3z+i) {— 131 ( 122+2 + I 12 + 1 ) + 132 ( 121+2 + Ill+l)} 

I Em I 

+ 7+—ry Ill +l(l3i+i) {—l3l(l22+2 + Il2+l) + l32(l21+2 + Ill+l)} • 

I Em I 
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On the other hand, the first two terms in (58) imply that = uj^-\-0{e) (ensuring 
consistency) and, moreover, that 


LJ 


fc+1 

1 

. fe +1 


/ 

Wi 

/ 

UJ2 






(60) 


which implies first order consistency in comparison to (44). Futhermore, = 

w*^ + 0(e) implies that the third and fourth terms in (58) cancel out at O(e^) order; 


thus, we realize that the relevant term in (58) at this order is just 

^2 / . .fc+1 


i: 


—CJ. 


,fc+l 




(I 3 . 


A‘+l 


(61) 


Using (60), we obtain 


1 


-a;r^(l3.ccr^) = - —(I3,a;f)^(l2icc2^ +Iiiccf) 
l-^m I 


1 


— — 7U!2{l3iUj^) {— 131 ( 122^2 +I 12 W 3 ) + 132 ( 121^2 +Iiia; 3 )} + 0 (e), 

l-^m I 


Wi"'‘^(l3iWy^) — --|^p-|(l3ia;f) ^(122^2 + 112 ^^ 1 ) 

+ jp—r^i (Isiwf) {— 131 ( 122^2 + I 12 W 3 ) + 132 ( 121^2 + IiiWj)} + 0{e). 
l-^m I 

Plugging these terms into ( |M] ), we obtain 

~ ItTUU {^ 22(l3i<^f )^(l2lW2 +I11W1) — Il2 )^ (122^2 +I12W3)} 

l-^m I 

— (1^^122 1^2 (l 3 iwf) {—131(122^2 +I12W3) +132(122^2 +I12W3)} 

|ilm I 

— ^1.,, (2^ ^12 (Isiwf) {—131(122^2 + 112^3) + 132(122^2 + I12W3)} , 

|llm I 

and 

^|J|- ( 2 ^ {l2l(l3iwf )^(l2lW2 +Ilia;3) — Ill (Isiwf) ^(122^2 +I 12 W 1 )} 

+ ^ J ^21 <^2 (Isi^i^) {~l3l(l22<^2 + Il2<^l) + I32(l21<^2 + Ill'll)} 

|llm I 

+ ( 2 ^ ^11 (Isiwf) {— 131 ( 122^2 + I 12 W 1 ) + 132 ( 121^2 + IiiWi)} . 

I Em I 


Now, substraction from the expressions of wi, 0)2 presented above, we obviously have 
|w(tfc + e) — ^ O(e^). Furthermore, the | factors in (58) prevent the O(e^) 

terms on the discrete and continuous sides to coincide. Regarding the multipliers, 
it is straightforward to see from ( |5^ that A/c+i — Wi (l 2 ia;f) — wf (IiiW^^) + 0{e), 
while X{tk + e) = A(a;^) + 0(e), where A(cli) is determined by (45). Hence, we see 
that A(t/j + e) — A^+i = Oq + 0{e) with 


On = 


|Ir„| 


((I32I21 — 131122)^2 + (I32I1I — l 3 lIl 2 )Wi) , 


which is different from zero, in general, and consequently the discrete multiplier is 
not consistent with the continuous one. This shows the first claim. 

□ 


Concerning the second, it is enough to apply theorem 5.2 
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Therefore, the variational integrators setting generates a second-order consistent 
method on the dynamical side (a fact which is interesting, since we are considering a 
first-order consistent discrete Lagrangian Id, and therefore we would expect a first- 
order consistent numerical method in the spirit of [33]) while it is not neccessarily 
consistent in the algebraic one. Needless to say, this is a drawback in the numerical 
scheme. However, due to the decoupling between the two parts mentioned above 
(we can obtain the values independently of the Afc+i’s), w e m ay perform, 

besides the reescaling, a discrete shift as described in remark 6.1 generating 


in consequence a (2, s) method for and therefore a (l,s) method for ([|). 

Unfortunately, such a shift cannot be understood in general as an alternate choice 
of the discretization map p, as long as the particular map is not just a rescaling 
anymore but involves nonlinear terms depending on the w variables. 

We illustrate these developments by the following discussion and plots. Before 
going into details, we point out that the energy is preserved in nonholonomic dy¬ 
namics along the solutions. This can be easily shown by direct calculation when 
we define the Lagrangian energy by — L{g, g), where we consider again 

the original setting expressed in equations 
we arrive at 


Taking the time derivative of 


dt dt 


dL 

W 


9 


dL 

dg^^ 


dL 


dg 


79 


dL 

dg^^ 


according to ([^. Thus, as long as g{t) is a solution of the Lagrange-d’Alembert 
equations, we have that p°‘{g)g = 0, and consequently El is conserved. In the 
reduced case, using similar arguments it is proven in [13] that the Euler-Poincare- 
Suslov equations preserve the reduced energy 




€£ 0 * 


where we consider the reduced Lagrangian l{^) ( jllj ). In our case, the reduced energy 
along the solutions reads U/(a)) = (with i = {1,2}). 

The preservation of this energy should be taken into account as a favorable property 
in the performance of nonholonomic integrators, as shown below. 

We consider an homogeneous rigid body with inertia matrix I = 
and initial values a;i(0) = 0.4 and W2(0) = 0.5 (with proper unities). We shall dis¬ 


1 

0.1 

0.2 

0.1 

1 

0.2 

0.2 

0.1 

1 

We shall dis- 


play the performance of an order (2,1) DREPS, in terms of (54), and the variational 
nonholonomic integrator (58) and (59), which is also order 2 in the dynamical vari¬ 


ables as proved in proposition 6.2 More concretely, the integrator DREPS corre¬ 


sponds to the midpoint rule in (54), namely Zi (w 

fc+i 


k ^k-\-l 


) = - (I3, 




and l 2 {uj^, ^ (Ig^ ^), where we recall that i = {1,2} 

(note as well that in this case there is no e dependence in Z 2 )- In order to achieve 
order 2 consistency with respect to (45), we determine A/^+i according to (54). 


On the other hand, the variational nonholonomic integrator (58) and (59) corre¬ 
sponds to the setup 0(w) = e/(w), (pdipj) = {a,Cj) and the retraction map r = cay. 
What we observe is that, for small time steps e, the behavior of both integrators is 
indistinguishable except with respect to the multipliers. As to be expected due to 
proposition |6.2] the nonholonomic variational integrator produces an inconsistent 
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discretization of the Lagrange multipliers. However, we display also the perfor¬ 
mance of both integrators over a big time interval, noticing that the variational 
integrator’s performance is much better in all respects, mainly with respect to the 
preservation of energy, where we observe a fast decay of the midpoint rule. More¬ 
over, we observe that the variational integrator even works quite well if the step 
size of the time discretization is taken to be relatively big. 
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Figure 1. In this figure we display the performance of the mid¬ 
point rule (DREPS(a;^, Afc+i; = 0, with inertia matrix I and 

initial values a;i(0) and W 2 ( 0 ) introduced above) for the nonholo- 
nomic rigid body with time step e = 10“^. The solid red line 
is obtained through a RK4 integrator (which we consider an ac- 
curete approximation of the continous nonlinear dynamics in the 
short term), while the blue dots represent the performance of the 
midpoint rule. The plots (a) and (b) correspond to the dynamical 
variables Wi, UJ 2 , while (c) displays the Lagrange multipliers A. On 
the other hand (d) shows the inconsistent multipliers generated 
by the nonholonomic variational integrator. Finally, (e) and (/) 
shows the preservation of the constraints and the energy Ei{Cj) up 
through round off errors, respectively. 
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Figure 2. This figure displays the comparison between the mid¬ 
point rule (the same as in Figures and the variational integrator 
(58), (59), for a time step e = 1^ = 1 (we recall that this inte¬ 
grator is also order 2 consistent in the dynamical variables). The 
former is represented by the green points and the latter by the 
blue ones, while the solid red line still represents the performance 
of a RK4. Variables wi (a), W 2 (5), A (c) and Ei (d) are displayed, 
while (e) represents the preservation of the constraints by the vari¬ 
ational integrator up through round off errors. We observe a better 
performance of the variational setup, mostly with respect to the 
preservation of energy, a fact which, considering the bigger time 
step, leads to the conclusion that its convergence to the actual 
solution is much faster and its long-term behavior will be more 
accurate. 
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7. Discretization of the Euler-Poincare-Suslov problem as 

PERTURBATION 


So far, we have focused on the discretization of the Euler-Poincare-Suslov equa¬ 
tions on 5'0(3), paying attention to its consistency order and D—preservation. In 
this section we return to any Lie group G (as long as it may be embedded in R") 
and its Lie algebra g, since the relevant results from above apply in this general 
case as well; in the last part we will particularize G = SO{3). 

As it is well-known, by discretizing the dynamics we introduce some discrepancies 
with respect to the continuous system even when the discretization is performed 
in some kind of structure-preserving fashion. It is needless to mention that this 
is a central and fundamental question for all kinds of numerical investigations, es¬ 
pecially concerning the long-term evolution of dynamical systems. Regarding this 
issue, we refer to m where a positive answer to the following question is given: Is 
it possible to embed a numerical scheme approximating the continuous-time flow 
of a set of autonomous ordinary differential equations (ODE) into the time evolu¬ 
tion corresponding to a non-autonomous perturbation of the original autonomous 
ODE? The positive result may be phrased as: Any p—th order discretization of an 
autonomous ODE can equivalently be viewed as the time—e period map of a suitable 
e—periodic non-autonomous perturbation of the original ODE (where e is the fixed 
discretization lenght); while the precise statement is: 


Theorem 7.1. Suposse that h S C"'(R”,R"), r > 1, and consider the autonomous 
ODE 


X = h{x). 


(62) 


Let F(t,x) be the fundamental solution of (62) satisfying F{0,x) = 0, and assume 
that there are an integer p> 1, a continuous function C : [0,cx)) —)■ [0, cx)) and a 
one-step difference approximation of step size e 


Xk+i = 4>{e,Xk), (0 < e < eo; fc e Z) 


which is consistent of order p, i.e. 

\(j)ie,x) - F{e,x)\ < G{\x\)eP~^\ 

Then, there exists a function d{e,t/e,x), as smooth as h and periodic int of period 
e, such that if G{t, s; e, a;)|^ G(s, s; e, x) = x, is the fundamental solution of the 
non-autonomous, e—periodic ODE 


X = h(x) e^d{e, t/e, x), 


(63) 


then 


G{e,0;e,x) = (l)ie,x), 


where G(e,0;e, •) : R" —>■ R" is the Poincare map (period map) for (63), corre¬ 
sponding to initial time s = 0. 

Following these ideas, in [20] the case of nonholonomic dynamics when Q = R” 
was explored, and the following theorem was proved. It is stated here for later use. 


Theorem 7.2. Let the matrix 


d^L 


be non-degenerate and positive-definite. 


dvidvl 

Then, any p—th order discretization of the nonholonomic ODE (evolving on D) 


^Do not confuse with the Lie group G. 
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may be embedded into the time evolution of a non-autonomous perturbation of the 
original nonholonomic DAE ([^ of the form 


g" = vl + ePdl{e,t/e,g,Vg), 


d^L 


• J _ 

dvldvi ® 


dL 

W 


+ ( 5 ) + )i(e, t/e, g,Vg), 


(64) 


= 0. 


Similarly, an analogous result can be obtained in the reduced setting. There are 
some differences though; first, we note that the Euler-Poincare-Suslov equations 
( [T7| ) are first-order rather than second, as are; second, the constraints (a“, C) = 0 
determine a subset C g of a linear space instead of a regular submanifold. The 
procedure to obtain the perturbation of the nonholonomic dynamics produced by 
a given discretization can be split into three steps, say: 

(1) Define an ODE evolving on g° from the Euler-Poincare-Suslov DAE p7| 
by projection, which we will call Euler-Poincare-Suslov ODE; 

(2) Choose appropriate coordinates in such equations to apply the theorem |7.1[ 

(3) Undo the projection process to recover the perturbed Euler-Poincare-Suslov 
DAE. 

Thus, we obtain: 


Theorem 7.3. Let 


+ ^aO-b + ^^db{e, tj e, ^), 


(65) 


be a non-degenerate and positive-definite matrix. Any 
p—th order discretization of the Euler-Poincare-Suslov ODE (evolving on ) may 
be embedded into the time evolution of a non-autonomous perturbation of the orig¬ 
inal nonholonomic Euler-Poincare-Suslov DAE 0 of the form 

91 

beZ Q^d 

( a “,0 = 0 , 

where we have used a local form of the equations where C®}, are the structure con¬ 
stants of g. 

Proof. Before going into details, we introduce some notation. We denote m = 
(niab) ■= while (to“^) = m~^ denotes its inverse (recall that (mab) is 

regular and positive-definite). With this, the first equation in ( [T7| may be rewritten 
as 

= fiO + Aam'^(a“)e, (66) 


where f^{f) := 


(1) Since fx°‘{g) = represent a set of independent one-forms spanning 

D° C T*G, it is easy to see that the a“ spanning (g°)° C g* are also linear 
independent. Therefore, we can decompose the algebra as 

g = g"" © (g^)-^, (67) 

with respect to the metric (mab)- Thus, from ( |66[ ) we can derive an ODE evolving 
on g° by eliminating the Lagramge multiplier. More particularly, if we take the 
time derivative of the nonholonomic constraints we obtain (a“,^) = 0 (note that 
a“ are constant), which, after replacing ^ by the right hand side of the equation 


(66), yields 


= ( 68 ) 
where (6“^) := ((a“,and (Cap) = It is possible to prove the in- 
vertibility of C using geometric arguments (see for instance [5D], lemma 2.4), or 
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just by arguing that m is full-rank and a“ constant-rank. Finally, we obtain the 
Euler-Poincare-Suslov ODE evolving on given by 

? = HO, (69) 

as long as (a“,0 = 0, where h{0 '■= f{0 + Aa(^)TO“^a“ and Aq(^) is defined in 

p). 


(2) Now, let us choose adapted coordinates with respect to the decomposition 
(67), say a = {a, a), where a = 1, ...,n, a = 1, ...,n — m (corresponding to g^) and 
a = n — m + 1,n (corrresponding to (g^ )■'■). Since g is a linear space, the choice 
of such adapted coordinates is trivial. Projecting (69) onto g^, we obtain 

r = HiO, 

where by ^ we are denoting O- Finally, we need to prove that two points belon¬ 
ging to g° can be connected up through a curve within g^ in one time-step e (a 
fact which is necessary for technical reasons in the proof of theorem 7.1). More 
concretely, given two points ^ 1,^2 G 0°, they can be locally written, with respect 
to the decomposition (67), as = (^“, 0) and 0 = iO, 6)- Consider the cut-off 
function 

Xo : El- —>■ [0,1] 

such that 

Xo(t) = 1 for t < 0, 

Xo(t) = 0 for r > 1, 

and assume that xo is real analytic for r ^ 0,1. Also, denote Xi(t) := 1 — xo('f). 
For instance, we could take 

Xo('f) = (1 + tanh (cot(7rT))) /2, 0 < r < 1. 


Defining 


Hit) ■= Xo{t/e) H + Xi(Ve) O, 

it is straightforward to show that c{t) = (^“(t),0) C g°, for t G [0,e], and that 
c(0) = and c(e) = ^ 2 - Now, we can apply theorem 1 7.1 1 to O = H{0, ensuring 
that any g—th order discretization can be viewed as the time—e map of a suitable 
e—periodic non-autonomous perturbation, namely 


e = hH0 + ePHe,t/e,0- 


(70) 


(3) “Unprojecting” (70) from g° to g we obtain the equation 

C =/i(0 + e^f^(eU/c 0 (71) 

as long as (a“,C) = 0- This basically finishes the proof, since this equation can be 
derived from a perturbed Euler-Pincare-Suslov DAE of the form 

i = fiO + Aa -k J(e, t/e, 0, 

(a“,O=0. 


(72) 


Roughly speaking, d{e,t/e,0) is the “component” ofd(e,t/e,^) along g°, as we show 
next. By a direct calculation, we see that the first equation in (72) leads to (©by 
eliminating the Lagrange multipliers, which gives 

Aa = XHO - m~^d), 

where Aq,(^) is defined in (68). Plugging this into (72) we finally obtain ©, 


including as well the identification d = m ^d— Cq ,/3 m ^ a'^ (a^, m ^d). 
This finishes the proof. 


□ 
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Of course, it would be ideal for a further analysis, if the perturbed Euler- 
Poincare-Suslov equations in (72) admitted a Lagrangian structure of some kind, i.e. 
if it resulted from a suitable reduction of the Lagrange-d’Alembert principle. The 
question, when this is true and how this is related to specific properties of both, the 
underlying discretization of the unperturbed problem as well as the interpolation 
procedure, is left as a subject of further study. 


Finally, when we restric to G = 50(3), it follows directly from theorems 5.2 
and|7.3|the next corollary. 


7.2 


Corollary 7.4. Under the regularity eonditions introduced in theorems \7.S\ and \7.^ 
ap—th order discretization of the nonholonomic ODE (evolving on D) and a p—th 
order discretization of the Euler-Poincare-Suslov ODE (evolving on so{S)^ ), related 
by the reconstruction equation = nfccay(ea)fc), obbey the following statements: 


( 1 ) 

( 2 ) 


the former may be embedded into a non-autonomous perturbation of the 


original nonholonomic Euler-Poincare-Suslov equations (651 , including the 
particular perturbation e^d, 

the latter may be embedded into a non-autonomous perturbation of the 


e”in(i-P)Jg and e™(bp) 


original nonholonomic DAE (64 1 , including the particular perturbations 


8. CONSLUSIONS 

We follow the ideas presented in [20j to study the discretization of the Euler- 
Poincare-Suslov equations on 50(3). First, concerning the consistency order of the 
discretizations corresponding to the unreduced and reduced settings, we find that 
when the discrete reconstruction equation is given by a Cayley retraction map both 
consistency orders are related to each other, remaining the same in the unreduced 
to reduced direction, and spoiled to first order (when the original discretization 
has order bigger than one) in the reduced to unreduced one. Furthermore, we study 
distribution preserving integrators, giving precise conditions under which this prop¬ 
erty is ensured from the discretization of the Lagrange-d’Alembert principle. On 
the other hand, we also show that this preservation may be shown for more general 
numerical schemes, and give a precise algorithm to achieve it based on the general 
reduced framework. These definitions and proved results are illustrated through the 
example of the Suslov problem on 50(3). We consider two different integrators, one 
of them based on the midpoint rule while the other is based on a variational scheme. 
As proved previously, both are order-2 consistent in the dynamical variables, but 
it is numerically shown that the variational one converges faster to the actual solu¬ 
tion and shows a better long-term behaviour. Finally, concerning the discretization 
understood as perturbation, we prove that any distribution preserving integrator of 
the Euler-Poincare-Suslov equations may be understood as a non-autonomous per¬ 
turbation of the continuous dynamics itself. In the case of 50(3) and considering 
a discrete reconstruction based on the Cayley map, the last result combines with 
the consistency bounds in order to provide the order of the perturbation produced 
in the unreduced dynamics by a reduced integrator, and conversely. 
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